AD-A210  481 


AR-005-377 


■■x~r 


SRL-C010-TM 


COPY  No. 


DEPARTMENT  OF  DEFENCE 


DEFENCE  SCIENCE  AND  TECHNOLOGY  ORGANISATION 

SALISBURY 

SURVEILLANCE  RESEARCH  LABORATORY 


SOUTH  AUSTRALIA  • 

f 

5 


TECHNICAL  MEMORANDUM 

SRL-001 0-TM 


SPECKLE  REDUCTION  IN  SYNTHETIC  APERTURE  RADAR  IMAGES 


i 

i 

l 


Technical  Memoranda  are  of  a  tentative  nature,  representing  the  views  of  the 
author(s),  and  do  not  necessarily  carry  the  authority  of  the  Laboratory. 


Approved  for  Public  Release 


Commonwealth  of  Australii 

V - '  MAY  1988 


89  7  25  046 


UNCLASSIFIED 


DEPARTMENT  OF  DEFENCE 

DEFENCE  SCIENCE  AND  TECHNOLOGY  ORGANISATION 
SURVEILLANCE  RESEARCH  LABORATORY 


AR-005-377 


TECHNICAL  MEMORANDUM 
SRL-OOIO-TM 


SPECKLE  REDUCTION  IN  SYNTHETIC  APERTURE  RADAR  IMAGES 
D.M.  McDonald 


,  SUMMARY 
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This  survey  reviews  noise  reduction  techniques  of 
particular  relevance  to  the  multiplicative  speckle  noise 
associated  with  synthetic  aperture  radar  images.  Where 
appropriate,  mention  is  made  of  the  Generalised  Synthetic 
Aperture  Radar  (GSAR)  software  package  developed  by 
MacDonald  Dettwiler  and  Associates  for  the  processing  of 
raw  SAR  data  in  a  standard  format  into  image  data. 
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1.  INTRODUCTION 

Synthetic  Aperture  Radar  (SAR)  has  been  developed  in  recent  years  to  improve 
the  along-track  resolution  of  sidelooking  airborne  radar.  Spaceborne 
platforms  such  as  Seasat  (1978)  and  the  shuttle-borne  radars  (SIR-A  and  SIR-B) 
have  allowed  the  monitoring  of  large  areas  of  the  earth's  surface  for  remote 
sensing  and  other  purposes.  Spaceborne  SAR  has  proven  to  be  complementary  to 
optical  and  infrared  imaging  systems  such  as  Landsat;  the  radar  system  is 
active  with  oblique  illumination  across  the  swath,  whereas  Landsat  is  passive 
with  a  constant  sun  elevation  across  the  image,  varying  with  the  season.  Thus 
only  the  radar  can  operate  day  and  night  and  at  any  latitude.  Radar 
wavelengths  can  penetrate  clouds,  allowing  imaging  of  normally  cloud-covered 
regions . 

Ease  of  interpretation  of  radar  and  optical  images  varies  with  the 
application;  for  example  the  oblique  illumination  of  SAR  enhances  faults, 
fractures  and  lineaments  (depending  on  the  look  direct icr.)  ,  whereas  st'-eams 
and  details  of  vegetation  tend  to  show  up  better  in  Landsat  images.  In  SAR 
ocean  studies,  internal  waves  are  detectable  through  their  modulation  of 
surface  waves,  with  the  potential  benefits  of  enhancing  our  understanding  of 
oceanic  processes  and  extracting  bathymetric  information. 

Because  of  the  coherent  nature  of  the  SAR  imaging  process,  detected  images 
suffer  from  the  presence  of  multiplicative  noise  or  speckle,  analogous  to  that 
observed  with  lasers.  This  speckle  contamination  confuses  the  interpretation 
and  classification  of  SAR  images  and  has  stimulated  much  research  into 
techniques  of  speckle  reduction.  Unlike  Landsat  multiband  images,  SAR  images 
to  date  have  been  singleband,  although  multiband  sensors  are  planned 
(eg  SIR-C).  Classification  of  radar  images  together  with  images  recorded  at 
different  times,  or  with  different  sensors  such  as  Landsat  have  been 
attempted . 

The  problems  of  segmentation  and  object  detection  are  closely  related  to 
classification.  Some  approaches  have  specifically  aimed  at  identifying 
geometrical  structures  such  as  edges  and  corners.  The  review  however  will 
concentrate  on  those  techniques  directed  at  improving  the  visual 
interpretability  of  SAR  images. 


2.  SPECKLE 


2.1  Introduction 

The  standard  model  for  the  image  restoration  problem  is  that  the  original 
intensity  distribution  f  is  acted  on  by  a  linear  system  with  point  spread 
function  h  to  produce  a  blurred  image  b.  The  recording  process  involves  a 
(possibly  non  linear)  transformation  S  and  introduces  additive  noise  n’ 
(usually  assumed  to  be  white  and  Gaussian). 

Thus  the  intensity  of  the  observed  image  g  is  given  by 


g  =  S(f--h)  +  n’ 


Many  techniques  are  available  for  restoring  images  contaminated  by  such 
additive  noise. 

The  first  important  characteristic  of  a  radar  speckle  model,  is  that  the 
original  intensity  distribution  f  is  corrupted  by  multiplicative  noise  n 
due  to  the  coherent  nature  of  the  imaging  process. 
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This  corrupted  distribution  takes  the  place  of  f  above,  and  is  acted  on  by 
recording  system  with  point  spread  function  h.  If  there  is  negligible 
additive  noise  and  the  system  is  linear,  then 


g  =  (f.n)-'-h. 


This  formula,  with  or  without  the  point  spread  function,  is  the  basis  of 
most  of  the  specific  speckle  noise  reduction  techniques.  The  point  spread 
function  leads  to  spatial  correlation  between  adjacent  pixels.  This  is  of 
particular  significance  when  the  final  image  is  oversampled,  for  example  to 
match  standard  map  scales. 

The  second  important  feature  cf  the  radar  speckle  problem  is  simply  the 
magnitude  of  the  variance  of  the  multiplicative  noise;  a  mean  of  1  with  a 
variance  equal  to  the  mean  squared  for  a  single  look  intensity  image. 
Improvements  in  interpretabiiitv  can  be  made  by  manipulating  the  look  up 
table  (LUT)  on  display,  (ie  the  table  relating  pixel  value  to  displayed 
brightness)  which  changes  the  effective  variance,  for  example  by 
implementing  a  square  root  function  (Matthews  et  al  1984) ,  or  a  high  pass 
filter  (Lim  and  Nawab,  1981).  The  combination  of  poor  SNR  and  the 
multiplicative  nature  of  the  noise  leads  to  particular  problems  in  edge 
detection  and  classification.  In  rhe  presence  of  speckle,  the  required 
spatial  frequency  bandwidth  to  recognise  even  a  binary  object  must  be 
several  times  the  spatial  frequency  bandwidth  of  the  object  itself 
(Arsenauit  and  April,  1986). 

The  terminology  of  image  processing  has  been  strongly  influenced  by 
photographic  recording.  It  is  convenient  to  refer  to  images  as  being 
either  amplitude  or  intensity  (amplitude  squared)  images.  Some  care  must 
be  taken  in  interpreting  the  literature  as  to  which  is  being  referred  in 
any  context.  Note  that  the  GSAR  package  produces  an  amplitude  image  as  its 
standard  product. 

2.2  Speckle  models 

Porcello  et  al  (1976)  modelled  a  typical  radar  scene  as  "a  collection  of 
specular  scatterers  superimposed  upon  a  background  of  diffuse  scattering 
surfaces".  The  diffuse  scatterers  form  the  speckle  background  and  the 
highly  coherent  microwave  transmitter  serves  as  a  source  analogous  to  a 
laser  in  a  laser-illuminated  optical  system. 

Goodman  (1976)  studied  the  speckle  observed  in  laser  systems  in  some 
detail.  He  viewed  speckle  in  statistical  terms  as  a  random-walk 
phenomenon.  He  considered  monochromatic,  polarised  radiation  u,  viz 


u (x , y , z ,  t )  =  A (x , y , z)  exp(i2irvt) 

where  v  is  frequency  and  A(x,y,z)  is  a  complex  phasor  amplitude 
ie  A(x,y,z)  =  | A(x ,y , z) | exp{  i8 (x ,y , z) } 


The  observable  is  the  irradiance  at  (x,y,z),  given  by 


I(x,y,z) 


t f iu^x-y.z;t)i2dt 


I  A(x,y ,z) | 2 


3 
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To  derive  the  probability  distribution  of  the  intensity,  an  analogy  was 
made  with  the  classical  random-walk,  problem. 

The  complex  amplitude  A(x,y,z)  is  regarded  as  being  made  up  of 
contributions  from  many  elementary  scattering  areas  on  the  rough  surface. 
Thus  the  phasor  amplitude  can  be  represented  by 


N 

A ( x , y , z )  -  ^  |ak|  exp  (i0R) 

k=l 


The  two  important  assumptions  which  are  now  made  are: 

( 1 J  the  amplitude  aR  and  phase  0R  of  the  contributing  elementary 

phasors  are  independent  of  each  other  and  of  the  amplitude  and  phase  of 
other  phasors. 

(2)  the  phases  of  the  elementary  contributions  are  equally  likely  to 
lie  anywhere  in  the  interval  (-ir.ir);  ie  the  surface  is  rough  compared 
with  the  wavelength. 

These  two  assumptions  allow  the  adoption  of  the  classical  random-walk  in 
a  plane  solution.  Provided  the  number  of  contributions  N  is  large,  then 
the  real  and  imaginary  parts  of  the  complex  field  at  (x,y,z)  are 
independent  zero-mean,  identically  distributed  Gaussian  random 
variables,  ie  the  speckle  can  be  described  as 


a(x,y )  =  aR(x,y)  +  ia^x.y). 

where  a_  and  aT  each  have  variance  a1 
R  I  a 

The  speckle  intensity  is  then 


I  =  I(x,y) 


a2R(x,y)  +  a2:(x,y) 


with  phase 


8  =  tan 


ax(x,y) 

aR  (x,y) 


The  intensity  or  irradiance  I  obeys  negative  exponential  statistics, 
ie  its  probability  density  function  is  of  the  form 


P(I) 


1 

Y  exp 


c-i/i) 


i  >  o 


where 
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I  =  2 a2  (Guenther  et  al,  1978) 

The  standard  deviation  of  this  function  is  equal  to  I. 

The  density  function  of  8  is  flat. 

Hence  the  probability  that  the  irradiance  exceeds  a  given  threshold  I 
is 


F(I  >  It)  =  exp  ( - I / I t )  ,  It  >  o 

Note  that  the  form  of  this  probability  density  function  (corresponding 
to  a  single  look  intensity  image)  explains  the  multiplicative  nature  of 
speckle  noise.  In  other  words,  for  a  uniform  scene,  the  standard 
deviation  of  the  intensity  equals  the  mean  intensity. 

If  the  amplitude  (rather  than  intensity)  of  the  image  is  considered, 
then  the  probability  density  function  becomes  a  Rayleigh  distribution. 
Small  intensity  values  will  be  spread  out  and  large  values  will  appear 
to  be  compressed  (if  the  result  is  scaled).  The  standard  deviation  (for 
one  look)  will  be  reduced  by  almost  a  factor  of  2  to  0.5227.  Converting 
from  an  intensity  display  to  an  amplitude  display  is  equivalent  to 
remapping  the  look  up  table  (LUT)  in  the  display  processor  as  previously 
mentioned.  Extensions  of  this  technique  will  be  discussed  later. 

The  contrast  ratio  is  a  convenient  measure  of  the  effect  of  speckle  on 
an  image.  It  is  given  by 


y  _  standard  deviation 
mean  value 

with  a  value  of  1  for  a  single  look  intensity  image.  It  is  the  basic 
aim  of  speckle  reduction  techniques  to  reduce  this  ratio,  for 
homogeneous  areas . 

The  approach  outlined  above  considered  an  array  of  many  scatterers 
randomly  distributed  in  a  cell  with  dimensions  large  compared  to  a 
wavelength.  The  central  limit  theorem  leads  to  a  Gaussian  probability 
distribution  for  the  total  field  arising  from  the  mutual  interference 
from  these  independent  scatterers.  Envelope  detection  leads  to  a 
Rayleigh-distributed  envelope  and  a  negative  exponential  intensity 
distribution.  The  n  th  normalised  moment  of  the  detected  intensity  I  is 
given  by  (Oliver,  1984) 


I 

n 


<In> 


<I>* 


n 


If  the  number  of  scatterers  is  reduced,  then  there  appear  additional 
terms  in  the  expressions  for  the  moments.  For  example,  Oliver  quotes 
the  second  normalised  moment  of  the  intensity  as  being 


Iz 


2d  -£)  + 


<ak> 

N<a2> 
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where  N  is  the  number  of  scatterers  and  a  the  scattering  amplitude. 

Various  models  have  been  proposed  for  the  scattering  process.  They 
incorporate  fluctuations  in  the  number  of  scatterers  or  the  scatterer 
cross  section.  The  independent  K-distribution ,  for  example,  arises  from 
a  negative  binomial  scatterer  density  distribution,  with  independent 
scatterers.  Oliver  proposed  a  model  using  a  T-lorentzian  cross  section 
fluctuation  (including  periodicity),  in  which  the  fluctuation  has 
r-distributed  statistics  and  lorentzian  spatial-frequency  distribution. 
This  model  includes  the  effects  of  correlation  between  scatterers  and 
the  finite  illumination  size,  and  reduces  in  the  appropriate  limit  to 
the  independent  K-distribution  model.  This  work  on  correlated 
K-distribution  models  for  representing  speckle  is  continuing  (Oliver, 
1985). 

Ouchi  et  al  (1987)  found  that  the  speckle  statistics  depend  on 
backscatter  cross  section  fluctuations  when  the  correlation  scale  of  the 
fluctuations  is  comparable  with,  or  exceeds,  the  SAR  resolution.  They 
found  that  in  addition  to  the  impulse  response  function  corresponding  to 
Gaussian  speckle,  the  autocorrelation  function  of  the  speckle  included 
the  convolution  of  the  autocorrelation  function  of  the  cross-section 
fluctuations  with  this  impulse  response,  referred  to  as  non-Gaussian 
speckle . 


3.  NOISE  REDUCTION  TECHNIQUES 

Noise  reduction  techniques  will  be  described  briefly,  together  with  the 
results  of  practical  comparisons.  These  comparisons  tend  to  be  subjective, 
when  applied  to  actual  images;  there  are  advantages  in  first  testing  the 
techniques  on  simulated  images.  However,  in  this  case  the  validity  of  the 
underlying  models  cannot  be  tested. 

3.1  Multilooking 

The  most  commonly  used  speckle  reduction  technique  is  that  of  multilooking. 
The  principle  behind  the  technique  is  that  the  sum  of  N  identically 
distributed,  real-valued,  uncorrelated  random  variables  has  a  mean  value 
which  is  N  times  the  mean  of  any  one  component  and  a  variance  which  is  N 
times  the  variance  of  one  component.  If  we  add  N  uncorrelated  speckle 
patterns  on  an  intensity  basis,  the  contrast  ratio  of  the  resultant  speckle 
pattern  is  reduced  by  JW. 

Uncorrelated  speckle  patterns  can  be  obtained  through  time,  angle, 
frequency  or  p>  larisation  diversity.  In  the  GSAR  processor,  azimuth 
multilooks  are  achieved  by  dividing  up  the  Doppler  frequency  domain  during 
azimuth  compression,  with  a  user-specified  maximum  overlap  of  successive 
looks.  Because  of  the  one  to  one  relation  between  along-azimuth  position 
and  Doppler  frequency,  this  is  equivalent  to  dividing  up  the  synthetic 
aperture  length. 

Lee  (1986)  has  compared  the  effects  of  using  amplitude  images  in 
multilooking  with  using  intensity  images  and  then  taking  the  square  root  to 
obtain  amplitude.  Lee  found  a  slight  advantage,  as  the  number  of  looks 
increases,  in  the  second  technique.  The  GSAR  package  in  fact  averages  in 
amplitude . 

As  shown  by  Porcello  et  al  (1976),  the  summation  of  N  uncorrelated  random 
variables,  each  with  the  identical  exponential  probability  distribution 

-i-  e  y/°o,  has  the  density  function 
o 
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P(U) 


1 

/  U  ^ 

N-l 

(ao/N)  (N- I ) ! 

iVN 

exp 

l".'* 

The  racio  of  che  mean  to 
as  pointed  out  earlier, 
of  the  amplitude  of  this 


the  standard  deviation  in  this  distribution  is  , 
On  taking  the  square  root,  the  standard  deviation 
speckle  noise  becomes  (Lee,'  1986) 


I  NT"  (N) 

Jr2(S+j) 


Note  that  for  a  single  look,  the  standard  deviation  is  reduced  from  1 
to  0.5227,  merely  by  taking  the  square  root.  This  is  almost  as  effective 
as  using  4  looks  in  intensity. 


There  are  some  limitations  with  the  multilooking  techniques.  The  most 
obvious  is  the  degradation  in  resolution:  a  synthetic  aperture  reduced  in 
length  by  a  factor  of  M  means  the  azimuth  resolution  is  degraded  by  the 
same  factor.  Further  problems  may  arise  if  there  are  specular  reflectors 
in  the  image,  or  other  strongly  angle-dependent  features.  Tomiyasu  (1983) 
modelled  the  speckle  in  a  pixel  in  terms  of  isotropic  scatterers .  He  also 
considered  the  effect  of  specular  or  quasi  specular  surfaces  comparable  to 
pixel  dimensions  or  larger,  where  the  response  is  highly  angle-dependent. 
Tomiyasu  (1984)  also  considered  the  effect  of  just  a  few  dominant 
scatterers.  He  found,  in  this  case,  some  correlation  between  the  spectral 
response  and  the  sub-pixel  texture.  Scivier  and  Orr  (1986)  have  developed 
an  adaptive  enhancement  technique  directed  at  the  situation  where  scene 
reflectivity  changes  during  the  integration  period.  Finally,  muitilooking 
must  be  incorporated  in  the  SAR  processor.  Most  of  the  techniques  to  be 
described  below  are  applied  to  the  processed  image. 

3.2  Image  domain  filters 

Li  et  al  (1983)  compared  multilooking  with  several  types  of  image  domain 
filter.  In  the  latter  case,  single-look,  full  resolution  images  are 
orocessed.  The  images  are  then  convolved  with  a  low  pass  filter  window 
which  eftectively  produces  a  weighted  average  of  several  adjacent  pixels. 
This  convolution  may  be  implemented  in  the  frequency  domain  before  final 
detection  (ie  on  the  Fourier  Transform  of  the  complex  image  pixels). 

Li  et  al  compared  five  filters,  designed  with  the  same  rectangular 
equivalent  widths 


REW  =  f  dx 

J  to  (o) 

The  filters  chosen  were 
(a)  Box  filter 

1  , 

0  ,  otherwise 


f,1*1 


,  -1<x<1 
,  otherwise 


(b)  Triangle  filter 


W(x)  = 
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(.  c  }  Sine  squared  filter  W'(x) 

(d)  Exponential  filter  W(x) 

(e)  Lorentzian  filter  W(x) 


The  test  was  performed  on  Seasat  scenes  and  an  equivalent  number  of  looks 
calculated . 


ES'L 


(EfP])2 

VAR [ P ) 


where  P  is  the  pixel  intensity. 

Improvements  were  found  using  this  measure,  the  degree  of  improvement 
depending  on  the  scene.  Different  filters  yielded  different  impulse 
responses  of  course,  but  the  question  of  the  impact  of  these  on  image 
interpretability  was  not  addressed. 

The  main  disadvantage  of  image  domain  filter  techniques  lies  in  the  amount 
of  computation  required  to  process  to  full  resolution  and  to  perform  the 
filtering.  A  further  1 imitation  lies  in  the  storage  requirements  for  full 
resolution  images. 

Guenther  et  al  (1978)  also  compared  several  techniques  including 
multilooks,  spatial  averaging  and  square  and  square  root  operations,  on 
simulated  digital  speckle  images  of  a  test  pattern.  The  effectiveness  was 
judged  by  a  set  of  observers  on  the  basis  of  minimum  detectable  contrast  as 
a  function  of  object  size.  They  found  no  change  for  the  two  non-linear 
filters,  but  however  a  subjective  improvement  in  noise  for  the  squaring 
filter.  Note  that  in  practice,  power  law  operations  for  visual 
presentation  can  be  implemented  by  modifying  the  look  up  table  in  the 
display  device,  and  are  thus  highly  efficient. 

Smith  (1978)  presented  a  noise  filtering  technique  which  in  principle  bears 
a  family  relation  to  the  IDF  techniques.  The  algorithm  consists  in  essence 
of  transforming  the  image,  applying  a  low  pass  filter  and  transforming  back 
to  the  image  domain.  Smith  employed  a  Hadamard-Wa lsh  transform  and  an 
adaptive  low  pass  filter  based  on  Students  t-distribution  to  remove 
statistically  insignificant  transform  spectra.  The  technique  was 
demonstrated  on  photographic  images. 

3.3  Additive  noise  reduction 

A  wide  range  of  linear  spatially  invariant  filters  have  been  used  in  image 
processing.  However,  these  have  been  developed  primarily  for  handling 
images  possibly  distorted  by  a  non-linear  but  spatially  invariant  process 
and  degraded  by  additive  noise  (eg  Chin  and  Yeh,  1983).  Such  filters 
include  weighted  and  unweighted  block  averaging,  Gaussian  and  various 
edge-detecting  filters.  The  smoothing  filters  tend  to  blur  edges,  and  all 
have  problems  when  applied  to  radar  images  degraded  by  multiplicative 
noise . 

For  the  simple  model  in  which  the  observed  image  comprises  an  ideal  image 
corrupted  by  multiplicative  noise  (ie  the  point  spread  function  and  any 
system  non-linearities  may  be  neglected)  then  homomorphic  filtering  may  be 
applied.  This  technique  preprocesses  the  observed  image  (in  this  case  by- 
taking  the  logarithm)  so  as  to  transform  non-additive  noise  into  additive 
noise.  After  applying  standard  noise  reduction  techniques,  the  inverse 
operation  is  applied  to  generate  the  estimate  of  the  ideal  image.  However 
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there  are  -'robiems  when  the  noise-degraded  image  is  operated  on  by  a  point 
spread  function.  Deconvolution  techniques  to  remove  this  effect  are 
available  but  their  performance  degrades  for  noisy  images  such  as  typical 
rt  images  (Frost  et  al  1982). 

The  median  filter  is  particularly  effective  in  removing  isolated  spot 
noise,  but  has  also  been  applied  to  images  corrupted  by  multiplicative 
noise.  Each  pixel  is  replaced  by  the  median  pixel  in  a  block  of  its 
neighbours.  This  non-linear  and  heuristic  filter  in  practice  preserves 
edges  much  better  than  the  linear  filters.  It  has  been  used  as  a  basis  of 
comparison  for  other  algorithms. 

Scollar  et  ai  (1984)  used  the  local  median  and  the  interquartile  distance 
instead  of  the  mean  and  variance  on  photographic  and  other  data,  for  spike 
removal  and  edge  preservation. 

Pratt  et  ai  (1985)  developed  a  pseudomedian  filter,  a  computationally 
efficient  approximation  to  the  standard  median  filter. 

Various  extensions  of  the  Kalman  filter  have  been  applied  to  images  by 
various  authors.  These  recursive  methods  deal  with  additive  noise,  and 
include  line-by-line  recursive  filters  as  well  as  extensions  to  two 
dimensions  (eg  Woods  and  Ingle,  1981).  These  models  use  a  semicausal  type 
of  model,  in  which  a  pixel  value  in  the  original  image  is  a  linear 
combination  of  previous  pixel  values  plus  an  adjustment  (input)  term.  The 
observed  image  consists  of  this  image  corrupted  by  additive  noise.  (Nahi 
and  Assefi  1972,  Beimond  et  al  1982,  Jain  1977,  Chen  1979.)  Lee  (1981) 
applied  the  Kalman  filter  to  a  linearised  approximation  to  the 
multiplicative  speckle  noise  case  to  derive  an  adaptive  local  statistics 
algorithm,  discussed  in  greater  detail  later. 

A  recent  publication  describes  an  adaptive  block  Kalman  filter 
(Azimi -Sadjadi ,  1987)  which  uses  an  autoregressive  image  model  with  an 
observation  model  that  includes  multiplicative  noise. 

3.4  Wiener  filtering 

Wiener  filters  have  been  applied  mainly  to  the  additive  noise  case.  They 
are  based  on  minimising  the  mean  square  error  (MMSE  filters)  between  the 
original  and  the  restored  image  and  require  some  knowledge  of  the  noise 
(and  true  signal)  characteristics.  Pratt  (1972)  presented  a  generalised 
Wiener  filtering  technique  in  which  the  filter  is  applied  to  a  transformed 
version  of  the  data  (for  example,  using  the  Fourier,  Hadamard  or 
Karhunen- Loeve  transformations).  Some  work  has  however  been  done  on  image 
restoration  in  the  presence  of  film  grain  noise. 

Walkup  and  Choens  (1974)  addressed  the  problem  of  additive,  signal 
modulated  noise,  ie 


r(x,y)  =  s (x , y )  +  cf ( s (x , y ) ] . n (x , y ) 


where  r(x,y)  is  the  noisy  image,  s(x,y)  is  the  ideal  image,  c  some 
constant,  f  is  usually  a  non  linear  function,  and  n(x,y)  is  a  zero  mean 
uncorrelated  noise  process. 


The  two  dimensional  Wiener  filter  has  a  spatial  frequency  domain  transfer 
function  of  the  form 
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W  (  u ,  v ) 


0  (U,V) 

rs 

0rr  (u,v) 


where  0  and  0  are  the  spectral  densities  of  r(x.y),  the  observed, 
rr  ss  y 

degraded  image,  and  s(x,y)  the  ideal  image  and  0  represents  their  cross 

spectral  density.  For  the  model  above  (with  zero  mean  uncorrelated  noise), 

it  turns  out  that  0  =  0  ,  and  the  resulting  filter  is  given  by 

rs  ss 

0ss(u’v) 

W(ll,V)  0  (u,v)  +  c2  0  (u  ,  V  )*0  ( U  ,  V  ) 

ss  s  s  nn 

i  i 


where  0  , 

s1sl  represents  the  spectral  density  of  s1(x,y)  =  f l s  <  x , y )  )  and  0 

is  the  spectral  density  of  the  noise  process.  For  white  noise  and 
stationary  image  statistics,  this  may  be  written  as 

*ss(u’V) 


W(u,v)  = 


0  (u,v)  +  c2Si  N 

ss  o 


where  Sj2  is  the  mean  squared  value  of  s^x.y) . 


This  filter  thus  has  a  higher  gain  at  spatial  frequencies  where  the  signal 
to  noise  ratio  is  high  than  it  is  where  the  ratio  is  low. 

Kondo  et  al  (1977)  extended  the  application  of  the  Wiener  filter  to  the 
case  where  the  image  is  degraded  by  a  system  point  spread  function  h,  ie 


I(x,y)  =  s(x,y)*h(x,y)  +  f [ s (x, y)*h(x,y) ] *n(x ,y ) 


=  s(x,y)*h(x,y)  +  s j (x,y) *n(x,y) 


Their  filter  is  of  the  form 


$  ( u ,  v ) .  H*  ( u ,  v ) 

ss  ' 

#  (u,v) |H(u,v) | 2  +  #  (u,v)*§  (u,v) 

O  O  h  J  11X1 

1  1 

where  H(u,v)  represents  the  Fourier  Transform  of  h(x,y).  The  relation 
between  this  filter  and  that  of  Walkup  and  Choens  is  clear,  when  h(x,y)  is 
replaced  by  a  6  function.  Kondo  also  considered  one  form  of  a 
multiplicative  noise  model: 


I (x,y)  =  n(x,y)» [s (x,y)*h(x,y) 


with  a  corresponding  filter 


#  (u,v) .H*(u,v)/n 

ss 


W(U’V)  #ss(u,v)|H(u,v)|2  +  ,4|2  4 


.  (u,v)*l#  (u,v) |H(u,v) j 

n|  n  n  ss 

i  x 


where  ni(x,y) 


n (x , y )  -  n 
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Note  that  this  form  of  a  multiplicative  model  is  different  from  that 
considered  in  coherent  radar  studies  by,  for  example  Frost  et  al,  viz 
I(x,y)  =  [s(x,y)*n(x,y]*h(x,y) .  As  will  be  discussed  later,  Frost’s  filter 
is  based  on  a  Wiener  approach. 

The  Wiener  filtering  work  described  above  has  been  directed  specifically  to 
photographic  studies,  in  which  the  multiplicative  noise  affects  the  signal 
after  degradation  by  the  point  spread  function. 

The  above  forms  of  filter  assume  that  the  image  statistics  are  spatially 
invariant.  This  assumption,  although  justifiable  for  photographic 
processes,  is  not  so  valid  when  applied  to  SAR  data.  One  extension  is  to 
assume  only  local  spatial  invariance,  updating  the  statistics  as  the  filter 
scans  across  the  image. 

3.5  Adaptive  filters 

3.5.1  Kansas  filter 

Frost  et  al  (1981,  1983)  considered  the  equation  as  discussed  previously 
I(x,y)  =  [s (x ,y)*n(x ,y) ]*h(x,y) 


where  I  is  the  observed  image,  s  the  (desired)  image  reflectivity,  n  the 
multiplicative  noise  and  h  the  system  point  spread  function.  -As  with 
the  Wiener  filter  approach,  the  MMSE  filter  for  stationary  image  data 
was  calculated  by  minimising  E((r(t)  -  1(c)  *  m(t))2].  This  approach 
leads  to  the  transfer  function 


M(f) 


nSr(f) 

Sr(f)*Sn(f) 


J, 

n 


1 

H*  ( f ) 


,  f  *  0 

,  f  =  0 


where  n  =  E{n(t)}  and  f  is  the  spatial  frequency.  Sr(f)  and  Sn(f)  are 
the  power  spectral  densities  of  signal  and  noise  respectively.  The 
filter  in  this  form  strictly  applies  to  a  homogeneous  area.  At  this 
point,  the  filter  is  recognisably  of  the  same  family  as  the  (Wiener) 
filters  derived  by  Kondo  et  al  and  Walkup  and  Choens . 


Frost  et  al  now  apply  a  standard  model  for  the  terrain  reflectivity, 
r(t),  regarding  it  as  an  autoregressive  process  with  an  autocorrelation 
R^ (t )  and  two  sided  power  spectral  density  of  the  form: 


Rrd) 


+ 


r 


2 


with 


S 

r 


(f) 


2o  2  a 
r 


a2  +  4ir2f2 


+  r  2  6  ( f ) 
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where  Che  values  of  the  parameters  a,  fl^2  and  r2  differ  for  different 
terrain  categories. 

The  model  for  the  multiplicative  white  noise  is 

R  (t)  =  o  2  6 (t )  +  n2 

n  n 

with 

S  (f)  =  a  2  +  n2  5(f) 

n  n 


Here  o^2  and  n2  are  scene  independent,  but  depend  on  the  sensor,  and  a 

varies  only  slowly  in  a  given  image  and  can  also  be  considered  constant. 
Substituting,  the  impulse  response  derived  is  given  by 


M'(t)  =  Kjcte'01^1 

where  a  depends  on  the  observed  image  characteristics  (o^/I)2. 
Frost  et  al  derived  a  direct  proportionality  of  the  form  a  =  K2(o^./I)2. 
Kj  and  K2  are  normalising  constants. 

This  filter  form  reflects  the  form  of  the  model  chosen  for  the  desired 
(ideal)  image  and  noise.  It  is  extended  to  cope  with  non  stationary 
statistics  by  adaptively  varying  a.  A  larger  value  of  a  implies  a 
narrower  impulse  response  and  less  averaging  than  a  smaller  value.  Thu; 
in  a  region  containing  an  edge  where  the  variance  o^2  is  large,  a  is 

also  large,  less  averaging  will  occur  and  the  edge  will  be  better 
preserved. 

The  filter  is  made  adaptive  by  updating  the  value  of  (a^/I)2  in  a  moving 
window  across  the  image. 

Frost  et  al  demonstrated  the  filter  on  a  Seasat  image  of  the  Goldstone 
Array.  They  concluded  the  filter  smoothed  noise  better  in  homogeneous 
areas  and  preserved  point  targets  and  edges  better  than  the  median 
filter . 

Freitag  et  al  (1983)  compared  the  Kansas  filter  on  agricultural  scenes 
with  the  unweighted  median  filter,  for  the  5x5  and  11  x  11  cases. 
Their  criteria  were  edge  detection  and  classification  accuracy.  There 
are  detail  differences  between  the  filter  implementation  in  the  Freitag 
paper  and  that  published  by  Frost  et  al.  However  the  principle  of 
narrower  filtering  where  the  image  is  noisy  is  maintained  in  order  to 
preserve  edges. 

Their  conclusions,  based  on  classifying  agricultural  fields,  were  that 
the  adaptive  filter  performed  better  for  high  contrast  edges,  but  that 
low  contrast  edges  were  blurred  by  averaging.  The  size  of  the  selected 
window  depended  upon  the  resolution  of  the  system  and  the  size  of  the 
homogeneous  elements  making  up  the  heterogeneous  scene. 

Panda  (1978)  developed  an  adaptive  filter  for  additive  noise  which  uses 
a  local  operator,  the  ’sample  window  autocorrelation'  to  identify 
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regions  which  are  noisy  and  do  not  contain  edges.  The  local  operator 
controls  an  heur istical ly  selected  filter  of  similar  operation  to  the 
Kansas  filter. 

i 

3.5.2  Lee  local  statistics  filter 

In  Jong-Sen  Lee's  adaptive  approach,  the  effect  of  the  point  spread 
function  is  neglected.  The  observed  image  is  then  expressed  as 


I(x,y)  =  s(x,y)*n(x,y) 


distribution 


the  multiplicative  noise  with 
(when  displayed  as  intensity). 


a  negative 


exponential 


The  observed  pixel  is  linearised  by  the  first  order  Taylor  expansion 
about  ( I ,  n) ,  viz 


I  =  ns  +  s(n-n)  +  n(s-s) 

The  filter  derived  from  minimising  the  mean  square  error  is  of  the  form 


ij 


+  K.  . 
ij 


(I 


ij 


n .  . 

ij 


s  .  .) 
ij 


where 


K.  . 
1J 


n.  .0.  . 

i)xij 


ij 


i? . 

ij 


'ij 


is  the  estimated  variance  of  the  noise-free  image,  given  by 


var (I .  .)  +  I . . 2 
ij  ij 


a  2  +  n?  . 
n  ij 


s  .  . 
ij 


2 


Lee's  filter  is  made  adaptive  through  estimating  the  'local  statistics': 
local  values  of  var(I)  and  I  are  estimated  from  the  image  (eg  a  window 
around  the  pixel  being  estimated)  and  o^2  is  estimated  either 

independently  (assuming  it  to  be  stationary)  or  by  some  appropriate 
algorithm  from  the  image.  The  filter  thus  amounts  to  a  linear  weighted 
sum  of  the  local  mean  and  the  image  itself. 


If  the  linearisation  mentioned  above  is  not  made, 
occurs  in  the  denominator  in  the  expression  for  . ; 


K.  . 
ij 


s‘.  o  ‘ 
ij  n 


+  n2 


n.  ,Q.  . 

ll_lj 


Q.  .  +  o2  Q. . 
ij  n  ij 


then  an  extra  term 
viz 


In  this  form,  the  filter  would  be  identical  to  that  derived  by  Kuan  et 
al  (1985)  for  multiplicative  noise.  The  approach  of  Kuan  et  al  may  be 
regarded  as  a  generalisation  of  Lee's  work,  and  can  be  applied  to  any 
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model  of  the  form  I  =  s+n'  ,  where  n'  is  unbiassed  and  possibly 
signal-dependenc  noise.  In  the  multiplicative  case, 


n'  =  (n- 1 ) *s 


and  thus  fits  their  model. 

The  advantage  of  Lee's  filter  is  that  is  does  not  require  a  model  of  the 
original  signal.  However,  artifacts  can  be  produced  (Lee,  1983b).  A 
disadvantage  lies  in  the  method  of  estimating  the  variances.  If  too 
small  an  area  is  used  to  estimate  the  observed  image  variance,  then  it 
may  happen  that  the  calculated  Qij  is  sometimes  negative.  This  can  be 
tested  and  set  to  zero.  In  addition,  if  the  noise  variance  is  spatially 
variant,  its  estimation  may  be  time  consuming.  Ideally  it  is  measured 
from  the  local  variance  of  a  flat  or  almost  flat  area. 

This  algorithm  does  not  smooth  in  areas  where  the  local  variances  are 
large  such  as  near  edges.  To  smooth  near  edges  while  preserving 
sharpness  Lee  (1981a)  extended  the  algorithm  essentially  by  redefining 
the  neighbourhood  for  determining  the  local  statistics,  according  to  the 
local  gradient. 

If  the  local  variance  exceeds  a  preset  threshold,  an  edge  is  assumed  to 
be  present  within  the  window.  The  window  is  divided  into  a  number  of 
possibly  overlapping  subareas  and  the  local  means  calculated.  A 
gradient  mask  is  then  applied  to  these  means  to  determine  the 
orientation.  The  area  used  to  calculate  the  local  mean  and  variance 
then  includes  only  those  subareas  on  the  centre  pixel  side  of  the  edge. 
The  increase  in  computing  time  is  kept  within  bounds  by  controlling  the 
threshold  above  which  the  extra  processing  occurs. 

Arsenault  and  Levesque  (1984)  combined  a  homomorphic  filter  with  the 
local  statistics  algorithm  (as  applied  to  additive  noise). 

Kim  and  Haralick  (1985)  compared  the  extended  local  statistics  filter 
with  weighted  median,  average  and  Gaussian  filters.  The  local 
statistics  algorithm  is  sensitive  to  the  estimated  noise  variance,  and 
Kim  and  Haralick  adopted  a  procedure  to  accommodate  the  observed  spatial 
variance  of  the  noise  of  their  images.  They  updated  the  noise  variance 
at  the  start  of  each  row  of  the  image.  The  technique  was  to  obtain,  for 
each  row,  n  local  variances  (one  for  each  of  the  n  pixels  in  the  row), 
order  them,  and  average  the  m  smallest.  Each  local  variance  is 
calculated  from  a  neighbourhood  around  the  centre  pixel. 

The  comparison  was  on  the  basis  of  noise  removal,  contrast  stretching, 
edge  enhancement  and  texture  preservation.  The  adaptive  filter  using 
local  gradient  and  local  statistics  removed  noise  well  without  contrast 
loss  or  edge  blurring,  but  tended  to  blur  images  in  homogeneous  regions, 
that  is  its  texture  preservation  was  only  'fair'.  The  weighted  median 
removed  noise  but  with  poorer  contrast  and  edge  blurring  properties,  but 
was  superior  in  preserving  texture. 

The  Gaussian  filter  removed  noise  without  contrast  loss,  but  smoothed 
edges.  These  three  filters  all  performed  better,  as  expected,  than  the 
average  filter. 

A  recent  publication  describes  an  adaptive  filter,  specifically  for 
unity  noise  variances,  based  on  writing  the  formula  for  multiplicative 
noise  as 
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I ( x , y )  =  s  ( x  ,  y )  +  (n(x,y)-l)»s(x,y) 


The  second  term  represents  an  unbiassed  signal  dependent  noise  term. 
The  filter  derived  is  of  the  form 


r.  .  *|i  - 


i . .  - 
2-J 


VAR(I  ), 


(I .  .  -  I  .  .  ) 
ij  ij 


The  filter  of  this  form,  using  the  local  statistics  calculated  over  a 
5x5  block,  was  compared  with  other  5x5  filters,  including  the  median, 
averaging  and  Lee's  adaptive  filter,  as  well  as  with  mult i looking .  This 
filter  was  considered  superior  on  the  basis  of  resolution  broadening  and 
the  derived  equivalent  number  of  looks.  The  filter  has  the  advantage  of 
comparative  computational  simplicity  (Nathan  and  Curlander,  1987). 
Extending  the  derivation  for  arbitrary  noise  variance  yield?  the 
extended  Lee  formula  (page  12). 


Mastin  (1985)  applied  5x5  filters  including  the  K-nearest  neighbour 
filter  of  Davis  and  Rosenfeld  (1978)  to  a  synthetic  image.  More 
recently,  Durand  et  al  (1987)  applied  ten  7x7  filters  to  a  SAR  image 
of  cropland;  they  preferred  the  extended  Lee  filter. 

Scivier  and  Orr  (1986)  adapted  Lee’s  approach  to  enhance  multi-look 
images,  aimed  at  the  case  where  scene  reflectivity  changes  during  the 
integration  period,  as  with  large  specular  targets  (Tomiyasu,  1983). 
The  essence  of  their  technique  is  to  obtain  local  statistics  information 
(the  local  mean  and  variance)  from  a  window  comprising  a  set  of  centre 
pixels  (one  per  look)  and  a  set  of  neighbouring  pixels  in  each  look. 
One  look  is  designated  the  'preferred  look'.  Thus  in  the  equation 


s..  =  I .  .  +  K . . (I . . 

1J  1J  1J  2-J 


V 


(Cf  P  11) 


I.  .  is  the  Dixel  value  in  the  preferred  look,  and 

ij 

I.,  (and  var(I..))  are  derived  from  the  mn2  pixels  in  the  nxn  window 
ij  ij 

over  m  looks. 


If  the  variation  between  looks  can  be  accounted  for  by  speckle  then  the 
enhanced  image  approaches  the  normal  multilook  image.  If  the  scene 
varies  significantly  between  looks,  the  result  approaches  the 
’preferred'  look. 

This  technique  was  demonstrated  on  a  Seasat  ocean  image.  The  results 
according  to  the  authors  suggested  a  greater  preservation  of  dynamic 
features  over  the  standard  multilook  method,  with  a  similar  degree  of 
speckle  reduction. 

3.5.3  Comparison  of  adaptive  approaches 

Frost  and  Lee  both  considered  a  multiplicative  noise  model,  and  Frost 
included  the  effect  of  the  point  spread  function  as  well.  Comparing 
Frost’s  approach  (neglecting  the  point  spread  function)  directly  with 
Lee’s,  then  similarities  and  differences  are  apparent. 
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Both  Frost  and  Lee  developed  a  filter  with  stationary  characteristics, 
and  then  make  the  filter  adaptive  by  using  local  image  statistics 
derived  from  a  window  about  the  pixel.  However,  their  development  of  a 
stationary  filter  differs.  Frost  essentially  uses  a  Wiener  filter 
approach  to  derive  the  equation 


M(u,v) 


n$  (u,v) 
s  s 


fss(u,v) 


He  then  invokes  models  for  the  noise-free  image  and  multiplicative 
noise;  the  noise  free  image  is  modelled  by  an  exponentially  decaying 
autocorrelation  function  for  the  signal  (plus  an  offset  to  account  for 
the  average  grey  level)  and  the  noise  is  modelled  as  white  uncorrelated 
noise  with  non-zero  mean.  This  specific  model  leads  to  a  particular 
form  of  the  filter  in  which,  by  manipulation,  the  observed  image 
statistics  feature  rather  than  those  of  the  noise-free  image.  This  is 
an  obvious  requirement  for  a  practical  filter,  but  the  cost  lies  in  the 
necessity  to  assume  a  specific  mathematically  tractable  image  formation 
model;  using  different  models  may  not  necessarily  lead  to  such  a 
convenient  final  form. 


Mote  that  Kondo  et  al  tested  their  Wiener  filter  (with  a  different  noise 
model)  by  using  a  real  image  to  which  they  added  the  effects  of  noise 
and  the  point  spread  distribution.  Thus  the  parameters  in  their  Wiener 
filter  equation  were  known.  As  would  be  expected,  they  found  superior 
performance  using  the  known  autocorrelation  function  of  their  original 
distribution  compared  to  assuming  it  to  be  white  uncorrelated  noise. 
They  did  not  try  models  such  as  that  used  by  Frost.  Walkup  and  Choens 
similarly  used  simulated  images  with  known  parameters. 

Lee  made  no  assumptions  about  the  noise  free  image  characteristics.  He 
approximated  the  expression  for  pixel  intensity  in  terms  of  an  expansion 
about  the  average  intensity, 


I .  . 
ij 


s..n..  +  n. .5 .  .  -  n..s.. 
iJ  ij  ij  ij  iJ  iJ 


and  applied  what  can  be  considered  a  Kalman  filter  approach  to  derive  an 
estimate  of  the  noise-free  pixel  value.  The  input  s^  is  operated  on  by 

a  linear  filter  with  transfer  function  n..  and  has  added  to  it  a 

ij 

zero-mean  noise  term  s..(n..-n..)  to  yield  I... 

iJ  ij  iJ  ij 

The  Kalman  filter  is  closely  related  to  the  Wiener  filter  (Sorensen, 
1970) . 

Comparing  Lee's  approach  with  Frost's,  the  final  form  of  the  filters 
derived  are  totally  different.  If  Frost's  model  is  valid  for  a 
particular  image,  then  the  use  made  of  this  prior  knowledge  (the  assumed 
autocorrelation  function)  in  principle  should  lead  to  'better'  results 
than  a  method  that  doesn't  use  prior  knowledge.  The  question  becomes 
one  of  the  validity  of  the  selected  model. 

These  two  adaptive  methods  share  similar  characteristics  in  that  the 
average  and  variance  in  a  local  area  are  used,  with  an  estimate,  derived 
in  some  other  step,  of  the  noise. 
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3.6  S igma  filter 

This  approach  is  based  on  the  sigma  probability  of  the  Gaussian 
distribution.  The  sigma  filter  selectively  averages  only  those  pixels 
within  a  given  xo  range  in  a  given  (2n+l)x(2m+l ;  window  while  excluding 

significantly  different  pixels,  with  a  modification  to  cope  with  spot 
noise.  If  a  range  is  specified,  then  for  the  multiplicative  noise 

case,  the  pixels  used  in  averaging  are  those  in  the  range 


I  .  .  -  2o  I . .  to  I . .  +  2o  I .  . 
ij  n  i]  ij  n  ij 

where  is  the  noise  standard  deviation.  The  unweighted  average  of  these 
pixels  replaces  the  centre  pixel. 

For  a  Gaussian  distribution,  95.5°o  of  random  samples  are  within  the  2 o 
range.  The  assumption  in  this  method  is  that  pixels  in  the  window  within 
this  range  are  from  the  same  distribution.  Thus,  substantial  edges  will 
tend  not  to  be  blurred. 

Isolated  spot  noise  may  not  be  removed  by  the  filter.  One  solution  is  to 
specify  a  threshold;  if  the  total  number  of  pixels  within  the  range  is  less 
than  the  threshold,  it  is  assumed  that  the  centre  pixel  is  an  isolated 
spot.  The  pixel  is  then  replaced  by  a  4-neighbour  average.  Alternatively, 
a  3x3  window  with  thresholds  can  be  used  in  another  filter  pass. 

.Vote  that  it  is  the  product  of  noise  standard  deviation  and  a  selectable 
parameter  that  specifies  the  range.  The  choice  of  this  parameter  is 
somewhat  subjective;  Lee  used  a  2o^  range  for  most  of  his  published 

examples,  with  several  passes  through  the  data  and  with  windows  from  3x3  to 
7x7.  Vote  that  for  specific  images,  where  details  of  interest  are  close  in 
grey  level  to  surrounding  pixels,  he  suggested  a  small  window  size  with  a 
reduced  range  (Lee  1983a,  1986). 

Lee  has  also  applied  a  biased  Sigma  filter,  in  which  pixels  in  the  upper 
intensity  range  are  averaged  separately  from  those  in  the  lower.  The 
centre  pixel  is  replaced  by  the  closer  average.  This  biased  filter 
enhances  contrast  and  sharpens  ramp  edges. 

The  Sigma  filter  is  a  relatively  efficient  technique.  Its  effectiveness 
compares  favourably  with  straight  averaging  and  median  filters,  although  a 
comparison  with  the  local  statistics  method  was  inconclusive  (Lee  1983a). 

Lee  (1983b)  also  compared  the  performance  of  the  sigma  filter  with  several 
other  filters  for  the  case  of  additive  noise.  These  included  the  gradient 
inverse  weighting  scheme  of  Wang  et  al  (1981),  found  to  be  inferior  in 
reducing  noise  variance.  This  filter  computes  the  absolute  differences 
between  the  centre  pixel  and  each  of  its  neighbours  in  a  3  x  3  cell.  The 
centre  pixel  is  replaced  by  the  average  of  its  own  value  and  the  value 
obtained  by  weighting  each  neighbour  by  its  inverse  absolute  difference. 
Thus  values  which  differ  greatly  from  the  centre  pixel  are  weighted  less, 
although  other  weighting  schemes  satisfying  this  principle  can  be  devised. 
The  principle  can  be  adapted  to  multiplicative  noise  by  modifying  the 
weighting  so  that  pixels  with  large  relative  differences  are  weighted  less. 
A  scheme  analogous  to  Wang  et  al's  could  replace  the  absolute  difference  by 
a  fractional  relative  difference  (less  than  or  equal  to  one)  based  on  the 
ratio  of  the  centre  pixel  to  the  other  pixels  (or  its  reciprocal). 
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Lee  also  considered  an  edge  preserving  filter  (Nagao  and  Matsuyama  (1979)) 
which  considers  a  5  x  5  neighbourhood  around  the  centre  pixel,  and  creates 
nine  overlapping  subregions,  comprising  eight  directional  regions  of  seven 
pixels  (including  the  centre  pixel)  and  one  3x3  subregion  around  the 
centre  pixel.  The  centre  pixel  is  replaced  by  the  mean  of  the  subregion 
having  minimum  variance.  Lee  points  out  that  the  sigma-filter  average 
could  be  used  in  place  of  the  mean.  However,  the  Nagao  filter  suffers  from 
shape  distortion  because  of  its  directional  subregions,  as  well  as  being 
more  computat ion- intens ive . 

3.7  Bayesian  techniques 

One  would  expect  that  prior  knowledge  about  the  original  image  and/or  the 
noise  distribution  should  lead  to  a  better  estimate  of  the  original  image 
than  methods  not  using  such  knowledge.  Representative  models  are 
available.  One  significant  group  of  reduction  techniques  are  those  classed 
as  Bayesian. 

Bayesian  methods  have  been  applied  mainly  to  additive  noise  images  (Hunt 
1977),  but  recently  Geman  and  Geman  (1984)  have  developed  a  technique 
applicable  to  the  multiplicative  noise  case  as  well. 

Their  technique  applies  to  degraded  images  described  by 


g  =  0 (H( f ) ) 


where  f  is  the  original,  g  the  observed  images.  H  is  the  point  spread 
function,  0  a  possibly  nonlinear  transformation,  N  an  independent  noise 
field  and  •  an  invertible  operation  such  as  addition  or  multiplication. 

If  the  above  equation  is  viewed  as  a  stochastic  problem,  then  one  must 
estimate  the  restored  image  given  only  the  recorded  image  and  some 
statistical  knowledge  of  the  noise.  The  posterior  conditional  density  is 
given  by  Bayes'  Law: 


P(f/g)  .  zL&m- mi 

p(g) 

The  MMSE  (minimum  mean  square  error)  estimates  are  the  mean  of  the 
posterior  density  p(f/g),  and  MAP  (maximum  a  posterior)  estimates  are  the 
mode  of  the  distribution.  ML  (maximum  likelihood)  estimates  correspond  to 
the  case  where  all  original  estimates  are  equally  likely. 

m2 

For  an  mxm  image  with  L  grey  levels,  the  number  of  possible  images  is  L  , 
which  rules  out  any  direct  search  for  the  optimum  restored  image. 

These  techniques  tend  to  require  iterative  runs  to  converge  on  an  optimum 
result . 

Habibi  (1972)  used  a  partial  difference  equation  to  realise  a 
two-dimens ional  recursive  model  (a  counterpart  to  the  one-dimensional 
Kalman  filter)  yielding  the  MMSE  estimate  for  additive  white  Gaussian 


noise . 
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Burkhardt  and  Schorb  (1982)  derived  an  MAP  estimate  using  the  Viterbi 
algorithm.  Naderi  and  Sawchuk  (1978)  developed  a  Bayesian  algorithm  based 
on  a  'noise  cheating  algorithm'  (Zweig  et  ai  1975),  which  they  showed  could 
be  justified  on  ML  grounds. 

The  technique  of  Geman  and  Geman's  MAP  method  is  applicable  to 
multiplicative-noise-degraded  images.  Their  noise  model  is  closer  to  Lee's 
than  to  Frost's,  in  that  they  do  not  consider  the  effect  of  the  point 
spread  function. 

Geman  and  Geman  adopted  a  stochastic  model  for  the  original  image, 
representing  it  as  a  Markov  random  field,  MRF  (equivalent  to  a  Gibbs 
distribution).  Pixel  grey  levels  and  the  presence  and  orientation  of  edges 
are  viewed  as  analogous  to  states  of  atoms  or  molecules  in  a  lattice- like 
physical  system.  The  posterior  distribution  is  also  an  MRF  with  similar 
structure  for  a  range  of  operations  including  blurring,  non-linear 
deformations  and  additive  and  multiplicative  noise. 

Gradual  temperature  reduction  in  the  physical  system  to  isolate  low  energy 
states  ('annealing')  corresponds  to  obtaining  the  MAP  (maximum  a 
posteriori)  estimate  of  the  image  given  the  degraded  image. 

Geman  and  Geman  developed  a  'stochastic  relaxation'  algorithm,  which 
generates  a  sequence  of  images  that  converge  to  the  MAP  estimate.  The 
sequence  evolves  by  local  changes  in  pixel  grey  levels  and  in  the  locations 
and  orientations  of  boundary  elements.  The  stochastic  relaxation  allows 
random  changes  that  decrease  the  posterior  distribution  as  well  as  changes 
that  increase  it.  This  reduces  the  problem  of  local  maxima,  which  can 
occur  with  deterministic  iterative- improvement  methods. 

For  a  neighbourhood  around  each  pixel,  a  set  of  cliques  is  defined  with 
associated  potentials.  These  potentials  are  defined  so  that  the  Gibbs 
measure  is  related  to  the  a  priori  probabilities  associated  with  the 
uncorrupted  image,  the  noise  process  and  the  sensor  characteristics. 
Calculating  or  measuring  the  clique  structures,  probabilities  and 
potentials  represent  a  significant  challenge  with  this  approach. 

Kuan  et  al  (1987)  developed  a  MAP  filter  using  local  statistics  (ie  a 
non-stationary  mean,  non-stationary  variance  or  NMNV  model).  Their  model 
requires  the  solution  of  a  cubic  equation  for  each  pixel  in  the  image. 

3.8  Geometric  filter 

Crimmins  (1985)  adopted  a  radically  different  approach  with  his  so-called 
geometric  filter.  This  algorithm  considers  the  gray  level  profile  of  an 
image  line  as  a  geometric  shape.  It  is  based  on  applying  a  single 
iteration  of  convex  hulling  algorithm  alternately  to  the  image  and  its 
complement.  It  is  essentially  a  one  dimensional  algorithm  applied  in  four 
different  directions  in  the  two-dimensional  image:  horizontal,  vertical 
and  the  two  diagonal  directions.  The  effect  of  the  filter  is  gradually  and 
iteratively  to  remove  narrow  valleys  and  towers;  it  tends  to  preserve 
spatial  information  . 

One  step  of  the  algorithm  considers  a  vertical  slice  of  the  image,  where 
the  height  represents  the  pixel  value.  One  line  thus  produces  a  two 
dimensional  graph,  in  which  pixels  in  the  umbra  (ie  below  the  image  curve) 
are  set  to  1,  and  the  rest  to  0.  Essentially,  in  one  step  of  the  convex 
hull,  a  pixel  value  is  changed  from  0  to  1  only  if  at  least  four 
neighbouring  and  contiguous  pixels  (out  of  eight  neighbours)  are  1. 
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The  algorithm  was  compared  advantageously  with  multi- looking  and  median 
filtering  on  synthetic  SAR  imagery.  It  is  used  at  ERIM  both  for  image 
presentation  as  well  as  preconditioning  for  further  computer  algorithms 
(Crimmins  1986)  . 


4 .  SUMMARY 

Speckle  smoothing  techniques  can  be  categorised  in  different  ways.  They 
include : 

(a)  general  methods  such  as  linear  filters,  primarily  applied  to  additive 
noise ; 

(b)  homomorphic  techniques,  in  which  a  multiplicative  degradation  can  be 

converted  into  an  additive  one,  allowing  the  application  of  the  general 

techniques  ; 

(c)  non-linear  filters  such  as  the  somewhat  heuristic  median  and  Sigma 

probability  filters,  based  on  a  window  around  the  pixel  of  interest; 

(d)  ~he  well-established  mult i- looking  techniques; 

(e)  image  domain  filters,  in  which  the  image  is  transformed  in  some  way,  a 
linear  filter  applied,  followed  by  transformation  back  to  the  image  domain; 

(f)  adaptive  techniques  based  on  the  image  statistics  in  a  window  around 
the  pixel  of  interest; 

(g)  Bayesian  techniques,  based  on  known  or  assumed  prior  knowledge  about 
the  original  image  and  noise  statistics;  and 

(h)  the  geometric  filter. 

The  standard  noise  model  takes  into  account  multiplicative  noise  only, 

although  the  refinement  of  including  the  point  spread  function  has  been 
made  in  the  Kansas  filter.  Prior  knowledge  in  the  form  of  known  or  assumed 
original  image  and  noise  statistics  has  been  included  in  some  approaches. 
Most  of  the  techniques  discussed  operate  on  the  processed  image;  however 
multi-looking  techniques  are  incorporated  in  the  GSAR  package. 

A  speckle  size  or  texture  can  be  observed  on  processed  images.  Some  recent 
work  has  considered  speckle  smoothing,  taking  advantage  of  this  spatial 
correlation  (Kuan  et  al  1987).  Such  an  approach  may  well  tend  to  be 

computation  intensive  but  represents  one  direction  of  future  work. 

The  multiplicative  nature  and  magnitude  of  speckle  noise  remains  a 
significant  problem  both  for  visual  interpretation  of  images  and  their 
classification. 
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